Quantum corrections to the dynamics of interacting bosons: beyond the truncated 

Wigner approximation. 
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We develop a consistent perturbation theory in quantum fluctuations around the classical evo- 
lution of a system of interacting bosons. The zero order approximation gives the classical Gross- 
Pitaevskii equations. In the next order we recover the truncated Wigner approximation, where the 
evolution is still classical but the initial conditions are distributed according to the Wigner transform 
of the initial density matrix. Further corrections can be characterized as quantum scattering events, 
which appear in the form of a nonlinear response of the observable to an infinitesimal displacement 
of the field along its classical evolution. At the end of the paper we give a few numerical examples 
to test the formalism. 



I. INTRODUCTION 

The huge interest to interacting bosons was stimu- 
lated by experimental advances in realization of ultra- 
cold Bose-Einstein condensates (BECs) 1,2,3,4 ' 5 . If atoms 
are a subject to an additional periodic potential com- 
ing from optical standing waves then they form perfect 
bosonic crystals with no dislocations, impurities or other 
defects. Moreover, the hopping amplitude between the 
adjacent cites of this crystal can be varied by orders of 
magnitude simply by tuning the intensity of the laser 
beam 3 ^. Even the sign and the strength of the interac- 
tion can be changed using an external magnetic field&t 
So one can experimentally address various phenomena 
without worrying about complications arising from the 
unwanted degrees of freedom. For example, in equilib- 
rium, the bosons can undergo a transition from a super- 
fluid to an insulator as the strength of the optical poten- 
tial is increased&SiiSiiiii^ii 3 .. This transition was directly 
observed in Ref. [j| . The simplicity of the system and rel- 
atively slow decoherence suggest a possibility to use the 
interacting bosons in quantum computingii 1 ^. A very 
appealing idea is to encode quantum information in dif- 
ferent sites of the optical lattice rather than in different 
internal atomic states 1 ^. In this way the measurement 
process and the manipulation of the information seem to 
be quite straightforward. 

Another huge advantage of the cold atomic systems is 
that one can address both theoretically and experimen- 
tally dynamic properties of the interacting atoms. In con- 
ventional many-body systems strong and non-adiabatic 
perturbations generically lead to fast damping due to 
strong coupling to various bath degrees of freedom. So 
if one is interested in quantum effects, it is extremely 
hard to work far from equilibrium. There are no such 
limitations, however, in the atomic systems. Indeed, in 
order to obtain extremely low temperatures where the 
atomic dc Broglie wavelength becomes comparable with 
inter-particle spacing, it is necessary to isolate the sys- 
tem to such extent that it can be really thought of as 
completely closed. One of the examples of the strongly 
out of equilibrium behavior occurs when the sign of inter- 



actions changes from repulsive to attractive, for example 
using Feshbach resonanceiii^iiS. The important issue 
for noncquilibrium problems from the theoretical point 
of view is how to take into account quantum fluctuations 
and go beyond classical or Gross-Pitaevskii (GP) equa- 
tions of motion. The latter can not adequately describe 
dynamical properties near instabilities simply because 
as the fluctuations grow in time, the initial wavepacket 
spreads very fast. The obvious generalization of the GP 
equations is the Bogoliubov's approximation^. The lat- 
ter, together with the perturbative treatment of inter- 
actions between quasiparticles, gives a consistent expan- 
sion of the partition function for equilibrium problems. 
On the other hand, in noncquilibrium systems it is not 
always possible to isolate a unique classical path. In par- 
ticular, if the classical motion can not be described by the 
linearized version of the GP equations, the whole Bogoli- 
ubov's expansion breaks down. In Ref. plj we considered 
a specific problem where neither GP nor Bogoliubov's 
theory would work. In that paper we studied a conden- 
sate, which was adiabatically driven to the point of the 
instability, and showed that because of quantum fluctu- 
ations the system evolves into a macroscopically entan- 
gled state. We argued that the correct dynamics can be 
reproduced within the truncated Wigner approximation 
(TWA) 22 i 23 i 24 i 25 i 26 i 27 . The idea behind TWA is that the 
evolution is still described by the classical GP equations 
of motion, but the initial conditions for the classical field 
are distributed according to the Wigner transform of the 
initial density matrix. We showed that the TWA natu- 
rally arises as the first quantum correction to the classical 
evolution and it gives asymptotically exact short time be- 
havior for any bosonic system. The main purpose of the 
present investigation is to show how one can go beyond 
the truncated Wigner approximation. 

There have been some works generalizing GP equa- 
tions by including the interaction of the condensate with 
excited bosons which are produced during the collapse of 
the condensate 28 29 . Another possible alternative to the 
GP and Bogoliubov's approximations is the conventional 
Keldysh technique 3 ^. As any diagrammatic expansion, 
this one relies on the smallness of interaction both for the 
classical and the quantum fields^. Therefore this tech- 
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nique is not suitable for strongly interacting systems near 
the classical limit. Also it is not suitable for the short- 
time dynamics, since the initial conditions are not explic- 
itly written, but rather absorbed into the quantum prop- 
agator in a complicated manner—. Quite different class 
of methods, which came from quantum optics, are based 
on the solution of the exact Fokker-Plank equation for the 
density matrix which is written in a coherent state rep- 
resentation 2 ^ 2 ^^. Because of the over-completeness of 
the coherent states, such a representation is not unique, 
neither is the resulting Fokker-Plank equation. For exam- 
ple, in the positive P-representation, the evolution of the 
density matrix can be mapped into the stochastic clas- 
sical dynamics 22-23 . Within the Wigner approximation, 
it is easy to reproduce the correct short-time dynamics 
and recover TWA, but going further in time becomes a 
tricky issue. Recently, there have been developed an ex- 
act stochastic representation of the evolution equations 
for the density matrix using the Fock basis, which might 
be preferable to the coherent state on o 33 i 34 . So far this 
theory is restricted to a particular class of two-body in- 
teractions and it is not completely clear (at least to this 
author) how to generalize it. Let us also point out that 
a semiclassical expansion of the dynamics for a system 
of interacting fermions has been discussed in a number 
of works by Filinov et. a/. 35 36-37 . There the authors 
suggested a perturbative scheme based on the integral 
representation of the evolution equations for the density 
matrix. Describing the dynamics in terms of the latter 
is certainly a valid concept. However, especially near the 
classical limit, it might be preferable to study the cor- 
rections to the evolution equations themselves. And this 
is exactly what we are going to do here. Our aim is to 
develop an expansion similar in spirit to that given in 
Refs. 35.36 37] . However, we will give a completely dif- 
ferent derivation, which is suitable for bosonic systems, 
and will give a simple intuitive interpretation of the ob- 
tained results. The other important difference between 
the two approaches is that we will not attempt to evolve 
the density matrix itself, all the dynamics will be ascribed 
to the observable operators. 

Another motivation for studying quantum corrections 
to the classical equations of motion is the possibility 
to take into account coupling to other external degrees 
of freedom, which are usually represented by a thermal 
bath. Though, as we mentioned above, the notion of a 
bath is not particularly useful for the condensates be- 
cause of their isolation from the environment, this cou- 
pling certainly provides a strong mechanism of decoher- 
ence in most of condensed matter systems. A standard 
way to add such a coupling into the picture is based on 
the representation of the bath by a set of harmonic oscil- 
lators satisfying certain physical properties. After inte- 
grating out the environment degrees of freedom one de- 
rives Langevin equations, which are essentially classical 
equations of motion with extra dissipative and random 
force terms 3 ^. Unfortunately, the Langevin equations can 
be rigorously derived only for noninteracting particles in 



a harmonic potential linearly coupled to the bath, i.e. in 
the limit where there are no corrections to the classical 
equations of motion due to quantum effects. So incor- 
porating the latter into the Caldeira-Leggett picture 3 ^ is 
another important theoretical challenge. Generalization 
of the derived evolution equations to open systems will 
be a subject of our future work. 

Before going into the actual calculations let us outline 
the main results of this paper. The quantum corrections 
manifest themselves in two ways: (i) the initial condi- 
tions for the classical equations of motion are distributed 
according to the Wigner transform of the initial density 
matrix; the classical observable is the Weyl symbol of 
the quantum operator, or equivalently its symmetrized 
version with fields substituted by c-numbers (see also 
discussion preceding |T5j) ), (ii) there are quantum scat- 
tering processes, which are represented as a nonlinear 
response of the observable to the infinitesimal trans- 
form of the fields along their classical trajectories. We 
would like to emphasize that (i) alone is equivalent to 
the truncated Wigner approximation. Amazingly the 
symmetrized quantum operators and the Wigner trans- 
form appear automatically in the path integral approach, 
where initially all the operators are normal-ordered and 
there is no obvious reason why the Wigner transform 
should ever emerge. For some situations it is sufficient 
to ignore the quantum scattering completely and to con- 
sider only the evolution along multiple classical paths. 
We would like to emphasize that any harmonic action (in 
particular the Bogoliubov's approximation) is completely 
described by TWA, i.e. by the classical dynamics with 
appropriate initial conditions. Moreover, to recover the 
Bogoliubov's approximation it is sufficient to linearize 
the classical GP equations around the stationary solu- 
tion. As we showed in Ref. JSJ, this linearization is not 
adequate for problems with unstable dynamics. Quan- 
tum scattering modifies the classical trajectories them- 
selves, but this modification is conceptually simple. In 
the first approximation there is only one scattering event 
(although distributed over the whole space and the time 
interval of the evolution), i.e., one has to add a small 
perturbation only once (for a local in time interaction) 
during the classical evolution and calculate the (nonlin- 
ear) response of the observable at the end. In the sec- 
ond order of perturbation theory there are two scattering 
events, etc. The whole picture thus resembles the per- 
turbative approach to the ordinary scattering problem 
in the Feynman path integral representation. Although 
here we mostly concentrate on the interacting bosons, 
the results are quite general. They are applicable to any 
type of evolution, where the classical equations of mo- 
tion arise as a saddle point of the action. Let us also 
point out that TWA is equivalent to taking into account 
all classical vertexes in the conventional Keldysh tech- 
nique 31 . Each quantum scattering event is equivalent to 
adding a quantum vertex. The perturbative approach 
we develop here should be recovered also in the Fokker- 
Planck master equation for the Wigner function^, see 
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also Refs. However, to obtain the nonlinear 

response of the observable, one must consider corrections 
to the distribution function in terms of functional deriva- 
tives with respect to the fields, so that the response is re- 
covered upon integrating by parts the product of the dis- 
tribution function and the observable. While this must 
be certainly possible to do, we will see that the functional 
integral derivation gives a very simple and elegant way 
of deriving the desired equations. 



II. QUANTUM CORRECTIONS TO CLASSICAL 
DYNAMICS 



Let us assume that a system of bosons is described 
by some hamiltonian: 7i(a,j, at, i), which in general de- 
pends on time and is expressed as a polynomial in cre- 
ation and annihilation operators and a. It does not 
matter whether a is a continuous in space or a discrete 
field sitting on a lattice. For a single band boson Hubbard 
model, which we keep in mind for specific illustrations, 
the hamiltonian reads: 



-J(o]oj+i +a] +1 0j) + — aja^a}^ - 1) 



the Plank's constant. To see the analogy explicitly let us 
define 



(3) 



Then, according to standard quantum-mechanical formu- 
las: 



(4) 



We had to scale the factor of vTV in © to have a unique 
limit for rij at N — > oo. It is easy to recognize in Q 
the conventional commutation relation between, say, the 
position and the momentum with l/N playing the role of 
h. The reason why the Plank's constant does not enter 
@ explicitly is that the phase 4>j does not have a classi- 
cal analogue, while the phase multiplied by H does have 
one; it can be related to the angular momentum in the 
rotating systems^. 

A time dependent expectation value of an arbitrary 
operator CI is given by: 



(!) 0(t)=^P( X )(x|TV T e^o«M^fi e - i /o«W^| x ) ( 5 ) 



where J is the tunneling constant between the neighbor- 
ing sites and U is the on-site interaction strength. Both 
U and J can explicitly depend on time. We also assume 
that the initial state of the system is given by a density 
matrix pq: 



Po 



= E p Wlx)(xl, 



(2) 



where \x) represents some basis and P(x) is the proba- 
bility to be in a particular state (this P would coincide 
with the Glauber P-functionp2S if {\x}} are the coherent 
states). If initially the system is in the pure state, the 
sum contains only a single term. For the Hamiltonian JQ, 
the classical limit is obtained as the number of bosons per 
site N tends to infinity. So we will develop the expan- 
sion in terms of l/N, where the latter plays the role of 



Here Tk t denotes time ordering along the Keldysh con- 
tour going from r = to r = t and then back to r = 0. 
The exponent of the operator is understood in the usual 
sense of an infinite product^: 



e ip H{T)dr = lim TT(i + i ArW(r,)), (6) 

q=0 

where r q = tq/Q and At = t/Q. The ordering Tk t 
requires that the multipliers corresponding to later times 
are placed closer to the operator CI. 

Note that contrary to the derivation of the evolution 
equations within the Keldysh technique^, we ascribe dy- 
namics to the operator CI rather then to the density ma- 
trix p. In the coherent state basis © reads: 



Cl(t) = / DafDa* f Da b Da% (a b o\po\afo)e- a h a f°+ a h a fi+ iH ( a /°> a fi) Ar ... e - a *fQ a fQ 
Cl(a*f Q,a b Q, t)e a f Q ab Q e~ a * b Q ab Q+a * b Q ab Q-i- m ( a >> q > a l q-i) At . . . e ~ a t,o a >>a_ 



(7) 



r 



Here 7i(a^b, Of b , r) and Cl(af : b, aj 6 , r) are the normal or- 
dered hamiltonian TL and the observable CI with opera- 
tors a and substituted by complex numbers and 
a*j b respectively. The expression above is intentionally 
written in the discrete form, since we want to take spe- 



cial care of the boundary effects. Instead of the fields 
a f and a b propagating forward and backward in time, it 
is convenient to introduce their classical {tp) and quan- 
tum (rj) combinations: ay = ip + rj/2, a b = i\) — r)/1. The 
names "classical" and "quantum" are not accidental since 
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in the classical evolution all trajectories are uniquely de- quantum fluctuations. After taking the limit At — > in 
hncd and the backward path should be exactly identical we derive: 
to the forward one. So a q (t) ^ comes entirely due to 



= J Dr]Dr]*DipDtp* - f |p|Vo + f ) (IW)* + ^^(*) " 

/ T Z^ g^ m j,*2n+i- m 2 2n m\(2n + l-m)\ 

n>l m=0 r r v ; 



Wo Q^m^2n+l-m 2 2n m\(2n + 1 - m) ! /' w 



where C(ip,i/j* , t) stands for the classical (GP) differen- 
tial operator acting on the field VK*)- 

tlWV1 4 t M . (9) 

Note again that the operator £ as well as the fields tp 
and ?7 contain spatial indices which we suppressed in (JSJ 
to simplify notations. The products like rf L in JBJ are 
understood as the appropriate sums: Equat- 
ing lEJ to zero gives the classical equations of motion. In 
particular, if TL is given by the normal ordered version 



of with aj and at substituted by tpj and tpl- respec- 
tively, then £j[ip, if)*, t] = is equivalent to the Gross- 
Pitaevskii equation: 

Let us explain in some detail how we arrive from Q 
to JHJ). There are several different contributions to the 
exponent in JJJ, which we call 5. The first one (Si) comes 
from the terms, which do not involve the hamiltonian TL: 



Q-i 



9 =1 



°/ 9 ( a /«+l~ a /<?) _ a bq( a bq- O-bq-l) 



+ a /o( a /i- a /o) - afeQ(afcQ- a&Q-i) - afa a w+ a/o(a6Q- (H) 



Here q denotes the discrete time while the spatial indices 
are suppressed. The first sum in the continuum limit 
transforms into the integral: 

Q-i 

^2 a /g( a /?+! ~ a ii) ~ o* bq {a bq -a bq -i) -> 

which under the substitutions a/ — * ip+r]/2, a b — ► iJj—rj/2 
and after integrating by parts becomes: 

jf K** (t) ^ - " (t) ^) + v/w#) - 

(13) 

In the continuum limit the first and the second terms 
after the sum in (|11J) clearly go to zero and the last two 



read: 

a *fQ( a bQ ~ a fQ) ~ a>l a b0 = -\ip \ 2 - M 2 / 4 

+1/2(VS% + - r(t)v(t) - \v(t)\ 2 /2. (14) 

Combining equations (flip - (|14J) we derive: 

-|^o| 2 -^-^ + i W o-^o). (15) 

The second contribution to the exponent of JJJ comes 
from the terms containing the hamiltonian. However, 
since all of them are proportional to At, the continuum 
limit becomes trivial and does not give extra boundary 
contributions: 
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S 2 = f dr(n(a f (r),a}(r),r) -n(a b (T),at(r),r)) = f dr 
Jo Jo 

WWr) + ^,^(r) + ^,r)-« W r)-^ ) ^)-2^ > r)) . (16) 
And finally the last step in our derivation is the expansion of the expression above in powers of r\: 

W(l&+2^+y)-W(^2^- T )=»7— ^— Hj ^ +^ ^ g^ 2 „ +1 _ m 22 » m , (2n+1 _ m) r ( 17 ) 



We intentionally separated terms linear in 77 and rf be- 
cause they enter the classical equations of motion, while 
the higher powers contribute quantum corrections. Com- 
bining (jnj, and (|T7|) we recover JSJ. 

Let us discuss how quantum fluctuations enter the clas- 
sical dynamics. From equation (JSJ we see there are two 
boundary and one bulk contributions. Thus, integrating 
out rj{t) results only in the modification of the observable 
operator: 

n cl = (Q(p + v y 2,^-^/2)). (is) 

Here the average is taken over 77 with the measure 
cxp(— \rj\ 2 /2). Note that the original operator fl entering 
must be written in the normal ordered form. It is 
easy to check that integrating out fluctuations according 
to (|18|) is equivalent to rewriting Q in the symmetrized 
form and substituting a and by ip and ip* . We also note 
that Sl c i coincides with the Weyl symbol of ft. Let us give 
a simple example illustrating this statement choosing 
to be a number operator: 

fl = a f a = -(a+a + aa f ) - ~. (19) 

According to (|18f) we have: 

n d = (p + rf/2, i, - n/2) = - \ (20) 

We see that f2 c ; obtained from the normal-ordered Q us- 
ing (|18|) is equivalent to substituting the symmetrized 
product of a and a) by ifnjj* in H19J1 • 

Another boundary contribution originates from the 
field 770 corresponding to the initial time. Because of the 
coupling to ipo, this fluctuation introduces a probability 
distribution for the classical initial conditions: 

p{ipo,tpo) = J dr)*dr]o{ipo ~ yHV'o + y) 

e -|-0o| 2 -3l»;o| 2 e §(»7o^a-»7oV'o) ) (21) 

Note that p(ipo,ipo) is nothing but a Wigner transform 
of the density matrix in the coherent state represent a- 
tiort 2 ^^ and therefore it is not a positively defined quan- 
tity. Sometimes p^Oi V'o) nas a we i r d nonlocal behavior 



and the semiclassical limit is achieved in somewhat non- 
intuitive way22, see also discussion given in Ref. [2]]]. If 
we ignore corrections to the classical equations of motion 
coming from the higher powers of the quantum field rj or 
the last exponent in @ , then the time dependence of the 
observable f2 will be given by: 

n(t)= /#S#oKV'o,V'o)^(t,^o,^), (22) 

where fi c ; (t, tpo, tpo) ^ s evaluated on the classical field ip(t) 
satisfying the initial condition ip(0) = i/jq. This expres- 
sion constitutes the so called truncated Wigner approxi- 
mation frequently used in quantum optics^. 

The deviation of the actual trajectories from the clas- 
sical paths comes from the bulk quantum fluctuations, 
which appear in the last multiplier in JSJ. Those can 
be interpreted as quantum scattering processes. Clearly 
they are nonzero only if there is an interaction between 
the bosons, i.e. H. contains terms of the order three or 
higher in the boson fields a and a'. For simplicity, let 
us restrict the following analysis only to the case of a 
two particle short range interaction, which is usually a 
good approximation for atomic gases (see 0). After the 
construction is clear, the generalization on other cases 
becomes straightforward. So 

Hint = ^jp a l a J ( a l a J ( 23 ) 
3 

We would like to point out that the unique classical limit 
at N — ► co is obtained when 

A(i) = NU(t) (24) 

is kept to be independent of TV— . Then the quantum 
scattering part of the action reads: 

S q = J* dr^- (r(r)v(r) + pr)rf{r)) \ V (r)\ 2 . (25) 

Because S q contains the third power (or in general also 
higher powers) of the quantum field 77, we can treat it 
pcrturbatively so that: 
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fi(t) = / #J#o2#o,V>o) 
, H) 2 > 

Q3 



W(r)5/,-(r) 



^ 3 



5 2 /^(r')a/ fe (r') 



where /j(r) is an infinitesimal shift of the field ipj(r): 

^(r) - ^(t) + /,(r), <#(r) -> ^(r) + /*(r). (27) 

Let us briefly outline the derivation of 12tjfl . Expanding 
in powers of r\ gives the following integrals: 



drf T d VT rf T m r£e 



(28) 



where H is the classical Hamiltonian of the nonlinear 
Schrodinger (or GP) equation, not to be confused with 
H: 



Next we use the change of variables: 



i>(Tq) = ip(r q -i) ~ iF(^(r,_i), r^fa-i), r,_i) + /(r g ) 

(30) 

and a simple identity: 



a; m e M2: dx = 2ir(-i) m 6 (m) (a). 



(31) 



Note that the transform from i/j(T q ) to f(r q ) is linear 
and gives no Jacobian for any interaction. The classical 
equations of motion are recovered from i|3U|l if / = 0. 
So nonzero / indeed corresponds to the deviation of tra- 
jectories from the classical ones. Clearly each term in 
the expansion in l|26l) gives an extra prefactor of 1/N 2 , 
which is the semiclassical parameter for the boson Hub- 
bard model. To see this, one has to keep in mind that all 
the fields must be rescaled as ipj — > VN ipj^, so that their 
expectation values become of the order of 1 and indepen- 
dent of N. Then each derivative with respect to / or /* 



jW 0V/(T)0/j(T) ") 

■ n c im,{f,n),i>*(t, {/=/*»,*) 
i 



(26) 



f,f*=o 



in (I2(j|) would bring an extra factor of 1/yN. The in- 
terpretation of H26f> and l|33[) below becomes transparent 
from the figure ^ The solid lines there denote classical 
trajectories and the cross represents a quantum scatter- 
ing event. The quantum corrections appear as a non- 
linear response to the infinitesimal displacement of the 
field, i.e. they reflect the rigidity of the classical motion. 
It is straightforward to generalize (|26J) to the interaction 
nonlocal in space or time. The only difference is that the 
field ip and its derivatives over / and /* would carry dif- 
ferent spatial or time indices. The variables / and /* in 



(29) 




FIG. 1: Schematic representation of the first quantum correc- 
tion to the classical evolution according to (12(i> and . Solid 
lines represent classical trajectories, the cross corresponds to 
a quantum scattering event. 



(|26|l are treated independently. In numerical evaluations 
it is more convenient to use: 







1 d 



i d 



d 



1 d 



i d 



df 2cW 2 93/' df* 2 cW 2 93/ 



Then l|26|l becomes: 



(32) 



fi(t) = / dVo^/wW-o^o) 



1 \ - 

3 



d 2 



d 2 



x 9^.(t) 



d 







N \d^ff(T) 93/|(r) 

fid W*, {/}),*) 



(33) 
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We would like to make a few comments about equa- 
tions and l|3^|) . In our specific example we used the 
number of bosons per well N as a dimensionless semi- 
classical parameter. As we discussed above, each term in 
the expansion in those equations brings an extra factor 
of 1/N 2 . The absence of H anywhere, may be confusing 
since the classical limit certainly corresponds to h — > 0. 
However, this should not be surprising since here and 
quite often in the atomic physics the Planck's constant is 
either completely absorbed into energies, which are mea- 
sured in Hz, or into time. We already argued above that 
the ultimate reason why h does not explicitly appear in 
our formulas is that the phase does not have a classical 
counterpart. If we use conventional observables like coor- 
dinate and momentum or angle and angular momentum, 
hr 1 will appear as a prefactor in the action justifying the 
saddle-point or classical approximation. In the same way 
the number of bosons per site N appears as a prefactor 
in the exponent of © after the rescaling ijj — ► y/~N ip and 
r\ — > \f~Nr\. So in general any expansion in the powers of 
r\ is in fact the expansion in powers of H. The other im- 
portant remark is that at small times the deviation from 
the classical dynamics due to the quantum scattering be- 
haves as 0(f(t)/N 2 ), where f(t) is some function of time, 
which vanishes at t — ► 0. This proves that the truncated 
Wigner approximation, where the quantum scattering is 
completely ignored gives the exact short-time asymptoti- 
cal behavior of the full quantum dynamics. This very re- 
markable result is to be contrasted with those obtained 
within Keldysh technique, where the short time scales 
are usually unaccessible. 



III. NUMERICAL EXAMPLES 

To illustrate the current approach let us consider some 
specific examples. Since the main purpose of this section 
is not to address particular problems, but rather numeri- 
cally test the formalism, we consider a simple case of two 
coupled condensates where it is straightforward to obtain 
the exact solution and thus to verify the accuracy of the 
expansion given in 1|26[) and (|33[l . Formally l|33|) . which 
we use in practice, is a multidimensional integral which is 
best evaluated using the Monte-Carlo methods. To find 
the first quantum correction we apply a small shift to a 
classical field at a random moment of time for each set of 
initial conditions, then follow the simultaneous evolution 
of the original and the shifted trajectories and at the end 
of the evolution calculate the response of the observable 
using finite differences. The first example we consider is 
related to the discussion given in our earlier worksSi*^S. 
Namely, we assume that the two condensates, described 
by the hamiltonian Q with j = L, R, were initially un- 
coupled with their wavefunction being a product of two 
number states: 



The sub-indices "L" and "R" correspond to the left and 
right sites respectively. It can be shownSi that the 
Wigner transform of the density matrix corresponding 
to the state (|34[1 is given by: 

KVo,</£) = 2e- 2 ^-l 2 +l^l 2 )i JV (4|^L| 2 )iiv(4|Vo fl J 2 ), 

(35) 

where Ljv(x) stands for the Laguerre's polynomial of 
the order N. Then suddenly at t = the tunneling 
was turned on and the fol lowing evolution of the sys- 
tem is studied. In Refs. |2ll40j | we showed that the 
Gross-Pitaveskii and the truncated Wigner approxima- 
tions are very good for sufficiently short time scales 
t <t c ~ N/X = J/U (we remind again that A = NJ/U is 
the dimensionless measure of the interaction and the clas- 
sical limit is achieved at N — + oo keeping A =const(iV)) 
and they break down completely for longer times t > t c . 
As an observable it is convenient to choose a scaled num- 
ber variance: 

n = W 2 (H aL ~ a « aj? ) ■ ( 36 ^ 

The classical counterpart of f2 can be found either from 
direct symmetrization of i|3t)|) or using (JTSJ): 

n =i = 4^(^-«H) a -8^2- ( 37 ) 

In equations (|lf)|> we can always rescale time t — > t/J 
so that the dynamics is completely described by a single 
dimensionless parameter A. For the illustration we will 
choose A = 1 corresponding to the intermediate interac- 
tion strength. In figure [21 we plot the resulting evolution 
of the number variance for the exact solution, truncated 
Wigner approximation and the first quantum correction. 
Although account of a single quantum scattering event 
does not considerably extend the domain of applicabil- 
ity of the classical description, it allows to determine the 
time t c , where the TWA breaks down without addressing 
the exact solution. Similarly, the second correction would 
show the time scale where the first one breaks down, etc. 
Relatively small extension of t c for this particular exam- 
ple should not be surprising. As we showed in Ref. |icj . 
the semiclassical description breaks down, because the 
phase difference accumulated between the different en- 
ergy levels becomes comparable to 7r and the discreteness 
of the spectrum becomes crucial. On the other hand in 
the semiclassical description the number and the phase 
are continuously distributed and it becomes very hard 
to approximate a discrete sum with essentially chaotic 
phases by a continuous integral. 

Another example where the discreteness of the spec- 
trum is crucial is quite opposite to that studied above. 
Assume that initially the system of bosons was in the 
nonintcracting supcrfluid state, which is characterized by 
the product of coherent states: 



10} = \N) L \N) R . 



(34) 



|o) = |Viv) lc |Viv) 2c 



(38) 
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FIG. 2: Scaled number variance defined in 1361 as a function 
of time (measured in the units of inverse tunneling) for the 
case of eight bosons per site. TWA includes only fluctuations 
at the boundaries of the time domain as described in the 
text. The first approximation also includes a single quantum 
scattering event. 



In this example we will not assume that the number 



of sites is equal to two, to avoid unnecessary complica- 
tions coming from the conservation of the total number 
of bosons. Then assume that the tunneling was com- 
pletely turned off. We will follow the evolution of the 
expectation value of a in a particular site. This example 
mimics the experiments by M. Greiner et. al. Q, where 
collapses and revivals of the condensate were observed in 
this way. In the absence of tunneling, the time evolution 
in a single site is described by the hamiltonian: 

H = A a t a ( a t a _ ly (39) 

By simple time rescaling we can always choose A = 1. It 
is easy to solve the corresponding Schrodinger equation 
in a number state so that 

(a(t)) = VNe N ^ U/N -^- u l N . (40) 

It is also straightforward to obtain analytic results within 
the truncated Wigner approximation and find the first 
few quantum corrections. Thus: 

a TWA (t) = VNe~^b (41) 

(i + w) 



ai (i) = a T wA(t) ^1 
a 2 (t) = a T wA(t) 1 



4/V2 

e 

47V 2 



it" 



it- 



127V2(i + ^_) 



it 6 

12N 3(l + JL) 
it 3 



7it b 



12iV3 (1 + JL) 247V 4 



2887V 4 (l + ^) 



7it 5 



96iV 5 (l 



240iV 4 (1 



2N ) 



192JV 6 (1 + ^)' 



240iV 5 (1 + iL) 



(42) 



(43) 



From the structure of 1411 - 143|) it is clear that the fur- 
ther quantum corrections correspond to the 1/iV Taylor 
expansion of l|4Uf) and this indeed can be explicitly veri- 
fied. However, both the truncated Wigner approximation 
(arwA(t)) and the quantum corrections fail to reproduce 
2-kN periodicity in time of the exact result I4H . So we 
can anticipate that the collapse of the condensate occur- 
ring at short time scales can be well described, while the 
revivals can not and this is indeed the case. The reason 
is that although the function (|4U|I is analytical in N, one 
needs to use a number of terms exponentially increasing 
with time to correctly reproduce (a(t)). On the physical 
level we can argue that the restoration of the coherence 
comes entirely from the discreteness of the spectrum, so 
that as in the first example, the continuous semiclassical 
expansion does not work at long time scales. The expec- 



tation value of the creation operator (a) as a function of 
time for the case of the four bosons per well is plotted in 
the figure 

Now let us consider few more examples, where the sug- 
gested expansion gives considerably better results. Imag- 
ine that initially two non-interacting coupled condensates 
are a subject to the oscillating in time interaction: 

X(t) = sin4i. (44) 

The frequency u — 4 is chosen to be exactly twice that 
of the Josephson-like oscillations in the noninteracting 
system, so that we can expect a parametric resonance. 
This process is thus equivalent to a resonant heating of 
the condensate. We would like to point out that the 
Gross-Pitaveskij approximation completely fails to de- 
scribe such a resonance, because the symmetric state is 
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FIG. 3: Expectation value of the creation operator as a func- 
tion of time (measured in units of 1/A) for the four bosons 
per well. The inset shows the same as the main graph, but for 
a longer time scale. The semiclassical expansion works very 
well for the short time dynamics, but fails to reproduce the 
restoration of coherence at a longer time scale. The second 
quantum correction I43H is not shown because it can not be 
distinguished by eye from the first one 1421 . 



a classical ground state for any strength of interaction 
A > — The number conserving Bogoliubov's approx- 
imation also would fail, because near the parametric res- 
onance we can expect (see figure 0} that a considerable 
fraction of particles will be excited. The resulting depen- 
dence of the relative number variance for eight bosons 
per well is plotted in figure Obviously, the truncated 
Wigner approximation gives a very good description of 
the evolution and the first quantum correction does even 
a better job so that it hardly deviates from the exact 
result. 

As the last example, we will choose an adiabatic evolu- 
tion of the ground state of the two coupled condensates 
with the increasing interaction: 



X(t) 



tanh St 
1-St' 



(45) 



Such a form of X(t) is chosen to mimic the tunneling 
amplitude J(t) exponentially decreasing with time, with 
the time being measured in the units of J (see Ref. 
Here 5 plays the role of the adiabaticity parameter. If we 
again take the number variance as an observable, then it 
will simply decrease with time. The question which re- 
mains however, is whether the truncated Wigner approx- 
imation will suffice to correctly describe the evolution or 
not. In the figure |3] we plot the number variance versus 
interaction for S = 0.05 and eight bosons per well. It is 
obvious that the TWA gives a systematic relative devi- 
ation from the exact solution, which grows in time and 
the first quantum correction considerably improves the 
agreement. 



FIG. 4: Relative number variance as a function of time (mea- 
sured in units of 1/J) for a resonant heating of initially nonin- 
teracting two coupled condensates with eight bosons per site. 
The interaction changes with time as X(t) = sin4i. 




1 2 3 
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FIG. 5: Relative number variance as a function of the di- 
mensionless interaction (A) for the interaction increasing with 
time according to I45H with 5 = 0.05. The truncated Wigner 
approximation gives a systematic deviation from the exact 
solution, while the first quantum correction makes the agree- 
ment with the latter nearly perfect. 



IV. SUMMARY 

Let us now summarize and discuss the derived re- 
sults. We developed a time-dependent perturbation the- 
ory around the classical evolution of the system of in- 
teracting bosons. We found the two types of correc- 
tions. The first one (which is equivalent to the trun- 
cated Wigner approximation) does not affect the evolu- 
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tion equations themselves but requires to take an aver- 
age over an ensemble of trajectories with initial condi- 
tions distributed according to the Wigner transform of 
the initial density matrix. The observable is the classical 
counterpart of the symmetrized quantum operator (or its 
Weyl symbol). Further corrections appear in the form of 
quantum scattering processes, which manifest themselves 
as a nonlinear response to the infinitesimal change of the 
fields along their classical evolution (see J2BJI). We would 
like to point out again that the widely used Bogoliubov's 
approximation, or more generally any expansion of the 
action up to the second order in the fields is entirely con- 
tained in the TWA. However, the latter is not limited to 
the processes, which involve only stable classical evolu- 
tion (see also discussion in Ref. |2ll|). 

Although we never discussed here the coupling to the 
thermal bath (or more generally to the external noise 
source), it is very straightforward to incorporate this into 
our picture. For example, in a simple representation of 
the bath as a set of harmonic oscillators, one obtains a 
dissipative term, which can be directly added to the GP 
equations and a quadratic term in the quantum fields, 
which is equivalent to a random force with a Gaussian 
distribution 38 . In general, there will be other terms in 
the action as well, but all of them can be treated pertur- 
batively in the same way as we explained irrespective of 
their origin. So they will result in some kinds of nonlin- 
ear responses to the classical, now stochastic, evolution. 
Note also, that the effects of the dissipation or random 
forces will wash out all the quantum scattering processes 
which occur a long time (longer then the relaxation time) 
prior to the observation. Clearly the response to an in- 
finitesimal perturbation will decay in time if we have such 
processes. So the theory can be extended to describe the 
evolution towards the steady states. One have to be cau- 
tious though in the low temperature limit, because gener- 
ically all the relaxation processes are frozen out and the 
time scale for the scattering leading to the equilibrium 
goes to infinity. This implies that as T — > more and 
more quantum scattering events have to be considered to 



correctly describe the equilibrium. So this theory would 
give some kind of a high temperature expansion. If one 
is interested in non-equilibrium dynamics, then the limit 
of the applicability of the perturbative expansion of the 
given order will be set either by the time of evolution 
or by the scattering time in a bath, whatever is shorter. 
Certainly further careful analysis of this approach for the 
systems interacting with bath is required and it is a sub- 
ject of the future work. 

Finally let us spend a few words on the numerical im- 
plementation of the calculations. The solution of the clas- 
sical equations themselves is very straightforward. Aver- 
aging over the initial conditions can be easily done with 
Monte-Carlo methods and the convergence time either 
slowly increases with the size of the system or saturates 
depending on a particular problem. But in any case it 
does not grow exponentially in size as would be the case 
for the full quantum solution. The quantum scatter- 
ing part requires evaluating a nonlinear response, which 
might be tricky for a nonlinear system of differential 
equations, but this part is also straightforward and well 
controlled numerically. We would like to point out that 
contrary to exact stochastic schemes, this method starts 
directly from the classical equations of motion, which 
may be very complicated themselves and shows how to 
add quantum corrections step by step. So it certainly 
must be very efficient if the quantum fluctuations are 
relatively weak. Thus the description of unstable "cat" 
dynamics considered in Ref. 21], which is very straight- 
forward in the present scheme can be hardly efficiently 
achieved using the stochastic equations. We believe this 
method is a competing alternative both to those based 
on the conventional diagrammatic technique, which is not 
very suitable for strongly interacting systems, and to the 
stochastic methods, which can usually deal with a limited 
number if degrees of freedom. 
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